A Linear Programming Approach to Routing Control in 
Networks of Constrained Nonlinear Positive Systems with 

Concave Flow Rates 


Heather Arneson a , Nicolas Dousse b , Cedric Langbort c 

a NASA Ames Research Center, Moffett Field, CA 94043-0001 

b Laboratory of Intelligent Systems, Ecole Polytechnique Federate de Lausanne (EPFL), 
EPFL-STI-IMT-LIS, Station 11, 1015 Lausanne, Switzerland 

c Department of Aerospace Engineering, University of Illinois at Urbana-Champaign, 
306 Talbot Lab/MC-236, 104 S. Wright St. Urbana, IL, 61801 


Abstract 

We consider control design for positive compartmental systems in which each compartment’s outflow rate is described by a 
concave function of the amount of material in the compartment. We address the problem of determining the routing of material 
between compartments to satisfy time-varying state constraints while ensuring that material reaches its intended destination 
over a finite time horizon. We give sufficient conditions for the existence of a time-varying state-dependent routing strategy 
which ensures that the closed-loop system satisfies basic network properties of positivity, conservation and interconnection 
while ensuring that capacity constraints are satisfied, when possible, or adjusted if a solution cannot be found. These conditions 
are formulated as a linear programming problem. Instances of this linear programming problem can be solved iteratively to 
generate a solution to the finite horizon routing problem. Results are given for the application of this control design method 
to an example problem. 
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1 Introduction 

Positive compartmental systems are popular models for 
describing interconnections of reservoirs whose dynam- 
ics are governed by conservation laws and natural pos- 
itivity and capacity constraints. Examples include au- 
tomobile or air traffic flows, job-balancing in computer 
clusters [4], or irrigation networks [3], to name just a few. 

A specific class of such compartmental systems, known 
as Eulerian models, has been used extensively in the air 
traffic management (ATM) literature [1, 2, 6, 7, 9, 10, 
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12, 13]. In an Eulerian model, the aggregate dynamics 
of groups of aircraft are modeled instead of the dynam- 
ics of each individual aircraft. As a result, the order of 
an Eulerian model depends only on the number of com- 
partments or sections used to describe the network of 
interest, but not on the total number of vehicles, which 
greatly reduces complexity in many cases. For more de- 
tails on Eulerian models and their comparison to La- 
grangian models, the reader is referred to [13]. 

Previous work in this area [1, 9, 10, 12, 13] has focused 
primarily on the use of linear models to describe the 
outflow of each compartment, or section, of the network. 
Such models describe the outflow rate of a section as 
depending linearly on the amount of material in that 
section. 

Here, we focus on concave outflow rate functions. Ex- 
amples motivating the use of this type of outflow model 
come from road traffic [8] and air traffic management 
research. The authors of [7] point out that, although 
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the outflow of a section of airspace will increase as the 
density of traffic increases, there is an upper bound on 
the outflow rate. At low density, flights are allowed to 
traverse a given section of airspace at their nominal 
speed and thus more aircraft in a given section results 
in a greater outflow rate from that section. At low traf- 
fic density, it is reasonable to assume a linear relation- 
ship between the number of aircraft in the section and 
the outflow rate of that section. However, aircraft must 
maintain a minimum separation for safety considera- 
tions, and thus, as traffic density increases, separation 
requirements cause the outflow rate to saturate. A non- 
linear, saturating outflow model is proposed in [7] and is 
shown to more accurately capture this saturating effect 
in dense traffic problems. In that work, the authors gen- 
erate an increasing, saturating outflow rate curve empir- 
ically through a discrete-event simulation which includes 
separation requirements. 

Our previous work [1, 2] has focused on variations of 
the problem of routing design for positive compartmen- 
tal systems to satisfy time- varying capacity constraints. 
The need to satisfy capacity constraints arises naturally 
in air traffic flow management problems. Airspace sec- 
tor capacity specifies the maximum number of aircraft 
that a trained human air traffic controller can safely 
route through the sector. Sector capacity depends on 
controller workload associated with traffic flow patterns 
in the sector. Additionally, the number of aircraft that 
can safely be routed through a sector of airspace at a 
given time can depend on the weather conditions in the 
sector at that time. Work has been done to estimate and 
predict sector capacities [11] in the presence of severe 
weather conditions. 

In [1] , we presented a solution to this problem for single 
destination networks with linear section outflow rates. 
The extension of this control design technique to the 
multiple destination problem is straightforward. Given 
linear outflow rates, a separate, decoupled network can 
be created for each destination. Routing solutions can 
be found individually for each sub-network using the 
technique presented in [1], 

We focused on single destination networks with nonlin- 
ear section outflow rates in [2]. In contrast to the lin- 
ear outflow problem of [1], extending this technique to 
the multiple destination problem is not straightforward 
because, with nonlinear outflow rates, sub-networks for 
each destination exhibit nonlinear coupling. This non- 
linear coupling makes the derivation of routing solutions 
more challenging. 

Here, we extend our former work and present a solution 
to the problem for a multiple destination network with 
nonlinear outflow rates. Rather than formulating con- 
straints, and subsequently a linear programming (LP) 
problem, to solve this problem directly, we make use of 
the solution for the multiple destination network with 


linear outflow rates. Additional constraints are imposed 
on the routing solution for the latter so that the result- 
ing closed-loop linear system behaves like the system 
with nonlinear section outflow rates. These constraints 
are nonlinear in control design variables and thus can- 
not be incorporated into the LP problem. To address 
this issue, we treat these constraint values as fixed and 
iteratively solve instances of the LP problem, adjusting 
the fixed values at each iteration. 

The multiple destination network with nonlinear outflow 
rates, control design objectives which ensure that basic 
network properties hold, and the formal problem state- 
ment are presented in Section 2. The derivation of the 
control design technique is presented in several stages in 
Section 3. First, the extension of the control design tech- 
nique of [1] for the multiple destination network with lin- 
ear section outflow rates, formulated as an LP problem, 
is described in Section 3.1. Control design for the multi- 
ple destination network with nonlinear outflow rates is 
then developed in Section 3.2. Constraints on the control 
input which force the multiple destination linear system 
to behave like the multiple destination nonlinear system 
are developed and incorporated as additional constraints 
in the LP problem of Section 3.1. A method to recover 
a routing strategy for the nonlinear system from a so- 
lution to this modified LP problem is given. Finally, an 
algorithm is presented to iteratively solve instances of 
this modified LP problem. In Section 4, an application 
example is given to demonstrate the proposed solution 
method. 

Notation: We are often concerned with vectors, scalars 
and elements of matrices which are indexed over some 
range. Therefore, for every positive integer N, we de- 
fine [ N ] = {1, 2, . . . , AT}. The cone of entry-wise non- 
negative vectors of dimension N is denoted by We 
write x > 0 to mean that x £ R+ , and x > 0 if it is in its 
interior. For all i £ [ N ], e* represents the i th canonical 
basis vector of . 

2 Problem Description 

2. 1 Model 

We consider positive conservative systems, which can be 
used to describe the flow of material through a network 
of interconnected reservoirs, or sections. 

In contrast with the models analyzed in [1, 9], we allow 
the outflow rate of each section to be a nonlinear func- 
tion of its state, which represents the amount of mate- 
rial present in the section. Nonlinear outflow rates are of 
particular interest in air traffic management problems. 
As discussed in [7], although the outflow of a section of 
airspace will increase as the density of traffic increases, 
there is a limit on the maximum outflow rate due to min- 
imum separation requirements between flights. Thus, we 
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are motivated by networks of air traffic which exhibit 
increasing, concave, saturating outflow rate functions. 
However, the control design technique presented here is 
applicable for concave outflow rate functions. 

The systems of interest describe the flow of material 
between interconnected “sections.” Sections may repre- 
sent, for example, reservoirs in an irrigation network, sec- 
tions of road in a traffic model, or volumes of air space in 
an air traffic flow management problem. Material flows 
between these sections making its way from one of sev- 
eral sources to a particular sink. A time- varying routing 
strategy is utilized to satisfy a particular performance 
objective. The performance objective of interest is pre- 
sented in Section 2.3. A control design technique is pre- 
sented in Section 3 which provides a method for gener- 
ating a routing strategy to satisfy this performance ob- 
jective. 


for alii £ [TV]. Note that, although physically relevant for 
some applications, the assumption that pi is saturating, 
that is 

lim fXi(a) < oo 

a— >■+ oo 

for all i £ [N ] , is not required for the application of the 
method and is thus not included in the above assump- 
tions. 

Assuming uniform distribution of material bound for 
each destination throughout each section, the portion of 
the outflow of section i bound for destination r is given 

by 

(2) 

Given the assumption that /q is differentiable at 0, and 
that Hi( 0) = 0, 


We begin by describing the dynamics of an /V-section 
network with R distinct destinations or sinks. In addi- 
tion to satisfying any performance objective, the routing 
solution must also ensure that all material reaches its in- 
tended destination. We aggregate material within each 
section of the network based on final destination. That 
is, assuming that there are R distinct destinations, we 
create R coupled sub-networks. Each sub-network de- 
scribes the flow of traffic through the IV-section network 
for a particular destination r £ [i?] . The state of section i 
of sub-network r at time t , which represents the amount 
of material in section i bound for destination r at time 
t , is denoted by x\{t). The total amount of material in 
section i at time t, denoted xflt) £ R+, is the sum of the 
states of section i for each destination, that is, 


lim 

Xj(i)— >■ 0 


fXj(Xi(t)) 

Xi(t) 


djjjj 

dxi 


(0) 


and thus the outflow rate given in (2) remains finite even 
as Xi(t) approaches 0. 

The fraction of the outflow of section i bound for desti- 
nation r routed to some subsequent section j at time t 
is specified by routing parameter /3£(t). To simplify no- 
tation, the set of routing parameters for sub-network r, 
namely { /3( ; (t)} t ., and the set of routing parameters for 

the full network {/^ (t)} . r , will be referred to as (3 r (t) 
and /3(f), respectively. 


R 

Mt) = '52 x i(t)- 

r= 1 

We define the state vector of each sub-network r £ [i?] as 


Some of the sections, referred to as “final sections,” are 
sinks through which material exits the network. The set 
of final sections for each destination r will be denoted by 
S 7 jr- It is assumed that the network consists of at least one 
sink for each destination, that is S r - F ^ 0 for all r £ [f?] . 


x r (t) = {x\(t), x r 2 {t), . . . , x r N (t)) T £ R+ . 

The state vector for the full network is 

x (*) = XXw 

r— 1 

= (xi(f),x 2 (t), . . . ,x N (t)) T £ M+. 

Motivated by our previous discussion, we restrict our- 
selves to outflow functions Pi : R+ —> K+ which satisfy 
the following assumptions 


w(0) = o, 

(la) 

f±i is differentiable at 0, 

(lb) 

( ii is concave, 

(lc) 


Network connectivity is specified for each destination. 
The subset of sections into which material in section i 
bound for destination r can flow is denoted by O'- . Since 
it may not be possible for material to reach a specific 
destination from any section in the network, we specify 
the set J\f r C [ N ] for all r £ [i?] as the set of sections 
i £ [N\ such that there exists at least one path from 
section j to a section in S'jr. That is, there exist sections 
i = ii,i 2 , - • - , i P such that £ 0” for l = 1,2, ... ,p — 
1, ii £ J\f r \Sjr for all l < p, and i p £ Sjr. We assume 
that the connectivity of all networks of interest satisfy 
0\ C A f r for all i £ [iV] , which ensures that material 
bound for destination r is routed into sections which 
are connected to destination r. Networks satisfying this 
property are said to be outflow connected. 

Material can always flow back into the section that it has 
just exited, therefore i £ OJ for all i £ J\f r and r £ [R]. 
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(c) Sub-network for destination 2. 



(d) Sub-network for destination 3. 


Fig. 1. Figures 1(a), 1(b), 1(c), and 1(d) depict the connec- 
tivity of the full network and sub-networks associated with 
destination 1, 2 and 3, respectively. The number in each box 
indicates the index of that section. Arrows indicate allow- 
able flow between sections of the network or sub-network. 
Recirculation is allowed for all sections. 

Physically, such recirculation corresponds to holding or 
slowing down the material moving through the section. 

Material is injected into the network by S sources, each 
one with flow rate d r s [t) for all s £ [S],r£ [R], and t > 0. 
The fraction of d r s ( t ) routed into section i is denoted by 
b r si , with 0 < b r si < 1 for all i £ M r , b r si = 0 for all 

i £ [lV]\Af r , and J2iLi Ki = 1 f° r all s £ [5], r £ [f?]. 
The inflow rate vector is defined as 

cf(t) = (cr 1 (t),d r 2 (t),...,d r s (t)f &r s , 

for each destination r £ [i?]. Inflow rates must satisfy 
d r s (t) > 0 for all s £ [5] and r £ [f?] . 

An example network consisting of 21 sections, three dis- 
tinct sinks and three sources is shown in Figure 1(a). 
Arrows indicate allowable flow between sources, sections 
and sinks. Recirculation is allowed in all sections. Fi- 
nal sections are 19, 20 and 21. The outflow of section 
19 flows into sink 1, the outflow of section 20 flows into 
sink 2 and the outflow of section 21 flows into sink 3 and 
thus = {19}, Sjr = {20}, and Sjr = {21}. Given that 
there are three distinct sinks in this problem, the full net- 
work is divided into three sub-networks. Sub-networks 
for destinations 1, 2 and 3 are given in Figures 1(b), 1(c) 
and 1(d), respectively. Notice that each sub-network has 
fewer than 21 sections. This is due to the fact that spe- 
cific destinations are unreachable from certain sections 
in the full network. For instance, in the sub-network as- 
sociated with destination 1, shown in Figure 1(b), it is 
not possible to travel from sections 18, 20 and 21 to sink 
1 given the section interconnection of the full network. 


These sections belong to the set [1V]\A/’ 1 and are omitted 
from the destination 1 sub-network. 

The dynamics of material in section i £ J\f r can be writ- 
ten as 

Xi(t) = V(x r ,x,P,fi,i) (3) 

for all r £ [i?] with initial state values 

a; r (0) = x r 0 , 

where x r 0 = x^ 0 , . . . , x r Nfi ) and 


T>(x r ,x,(3,n,i) = < 


+ (t)) 

+ £f=i^<m if * £ RT r , 


0, otherwise. 


The initial state must satisfy (xq); > 0 for all i £ J\f r 
and (xg)i = 0 for all i £ [_/V]\A/’ r . 


We denote the solution to (3) with initial state a; r (0) = 
x r 0 for all r £ [R] and inflow rates d r s ( t ) for all s £ [S'] and 
r £ [R], under the particular choice of routing strategy 
/3(t), as x /3 (t). 


2.2 Basic Control Design Objectives 

The model introduced in Section 2.1 describes the flow 
of material through a network of sections. In order to 
be physically meaningful, this model must satisfy the 
following constraints: 

Positivity: Given initial state Xq > 0 for all r £ [R] 
and inflow rate d r (t) > 0 for all r £ [1?] and t > 0, the 
resulting state vector satisfies 

x r (t) > 0 

for all r £ [R] and t > 0. 

Conservation: Constraints C(f3 r (t)) are defined for all 
r £ [i?], t > 0 by 


Vm,£AT (4) 

Pij(t) = 0, Vj £ N r \O r i ,Vi £ Af r (5) 

£ W) = i, v; e .Arvs> (6) 

jeoi 

£/%•(*)<!, Vi €S£. (7) 

jeoi 
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Note that constraints C(f3 r (t )) do not include terms 
/3 ij(t) for i £ Af r or j ^ Af r . It is not necessary to en- 
force constraints on these terms as they do not appear 
in the dynamic model (3). When required, individual 
constraints included in the set C(/3 r (t)) will be referred 
to using subscripts. For example, constraint (4) will be 
referred to by Ci(/3 r (t)). 


when Xj(i) = 0 given assumptions (la) and (lb). Thus 
x p (i) > 0 whenever x p (t) = 0 and we can conclude that 
x r (t) > 0 for all r £ [1?] and t > 0 and consequently 
x(t) > 0 for all t > 0. □ 


2.3 Performance Control Design Objective 


Positivity ensures that the state of each section of ev- 
ery sub-network, which represents a physical quantity, 
remains non-negative at all times. Physically, the Con- 
servation requirement ensures that material leaving a 
given non-final section enters another section in the net- 
work and is thus conserved. When designing the control 
input, care must be taken to ensure that the system sat- 
isfies these two conditions. 

Conservation is satisfied when C(/3 r (t)) hold for all 
r £ [A]. Given non-negative routing parameters, that 
is, routing parameters satisfying Ci(/3 r (f)), constraints 
C 3 (/3 r (f)), ensure that material exiting a non- final sec- 
tion is routed to a subsequent section in the network. 
Note that constraints C±{f3 r (t)), defined for final sec- 
tions only, allows material to exit the network. 

Again, given non-negative routing parameters, the Pos- 
itivity constraint is satisfied when the section outflow, 
/ii for all i, satisfies certain conditions, detailed in the 
following claim. 

Claim 1 Let initial condition Xq > 0 for all r £ [A], 
inflow rates d r s {t) > 0 for all s £ [S], r £ [fi], t > 0, 
and routing parameters ftljft) > 0 be given for all i,j £ 
Af r , r £ [1?], t > 0. Then , if each outflow rate function p,i 
satisfies assumptions (la) and (lb) and 

Pi(a) > 0 for all a > 0, 

for alii £ [IV], the solution of system (3) satisfies x r (t) > 
0 for all r £ [A] and t > 0 thus Positivity holds, and 
consequently x(t) > 0 for allt > 0. 


The objective is to design routing parameters (3(t) to en- 
sure that system (3) is Positive and Conservative , and 
the solution x l3 (t) satisfies the specified capacity con- 
straints. We assume that section capacity updates would 
be issued at regular intervals. Thus, we assume a piece- 
wise constant capacity constraint, as proposed, for ex- 
ample, in [5]. 

Given initial conditions and capacity constraints, a rout- 
ing solution may not exist or, it may not be possible 
to find a feasible routing solution using the proposed 
control design technique. If this is the case, rather than 
not return a solution, we would like to return a solution 
which minimizes, according to some metric, the viola- 
tion of the given constraint. 

The problem of interest can formally be stated as follows: 

Problem 1 Given a piecewise constant vector-valued 
capacity constraint profile c (t) on a finite time horizon 
[0,1V], find an adjusted capacity constraint c(t) > c (t) 
and routing strategy (3(t) ensuring that system (3) sat- 
isfies Positivity and Conservation and the capacity 
constraint condition xPft) < c (t), is satisfied for all 

t £ [0, H] while minimizing fj 1 (c (t) — c (t)) dt. 

Here, we assume that network inflow rates and inflow 
routing are fixed. These values could be incorporated as 
control design variables, in which case additional con- 
straints and objective terms would be required to force 
material through the network. 

3 Control Design 


PROOF. Given that Xq > 0 for all r £ [A], in order 
to show that x r {t ) remains positive for all r £ [A] and 
t > 0, it is sufficient to show that for any i £ Af p and 
pe[R] 


*?(*) = o 

x r (t) > 0,Vr £ [A] 


=> &i(t) > 0 . 


With xf(t) given by (3) and assuming that Xj(t) > 0 
for all j £ A C and r £ [A], and, in turn, x;(t) > 0 for 
all l £ [IV], it is clear that the only potentially negative 

component of xf(t ) is — /q (xfit)). However, when 
%?(t) = 0, nfixfit)) = 0, since is finite, even 


Our previous work [1] focused on the solution of this 
problem for a single destination network in which the 
outflow rate of each section depends linearly on the 
amount of material in that section. This solution method 
can easily be extended for application to networks with 
multiple destinations and linear outflow rates. The con- 
trol design technique presented here is based on this 
solution to the multiple destination, linear outflow rate 
problem. The significant advancement of this technique 
is the derivation of additional constraints and a method 
for generating state-dependent routing parameters for 
the nonlinear system based on the routing solution for 
the linear system. Together, these constraints and the 
routing parameter generation method ensure that the 
resulting routing parameters satisfy the Conservation 
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constraints and that the state of the closed-loop nonlin- 
ear system satisfies the capacity constraints. 

In Section 3.1, the control design method of [1] is ex- 
tended to the multiple destination problem with linear 
outflow rates. This is a relatively simple and straightfor- 
ward extension, and is presented only as an intermediate 
step in the solution to the multiple destination nonlin- 
ear problem. The proposed solution to Problem 1 is then 
presented in Section 3.2. 

To begin, we construct a continuous capacity bound c (t) 
in a piecewise manner using a subdivision of the time 
horizon of interest, [0,77], into K intervals of size At 
(chosen to be a divisor of 77), where K = It is as- 
sumed that the capacity constraints c(t) are constant 
over intervals of length At. The control design prob- 
lem then becomes one of finding a routing strategy /3(t) 
which ensures that x /3 (t) < c(t) where c (t) < c (t), when 
possible. If this is not possible, we would like to choose an 
adjusted capacity constraint c(t) such that c(t) < c(t), 
c(t) > c(t) and f Q (c (t) — c(t)) dt is minimized. 

We will be dealing with discontinuous functions, there- 
fore for any function g we define 

ffOfc) = t 1 ™ 9{t) and g(t~) = ^lim g{t). 

t>tk t< t k k 


parameters y are analogous to x r , x, and /3, defined in 
Section 2.1. 

The outflow rate of section i is fflyfli)) = Aifi where 
constant r, > 0 is the average traversal time of section 
i. The dynamics of section i with destination r is 

vUt) = ' D (y r :y^v,f,i)- (8) 

Assuming a uniform distribution of material bound for 
each destination throughout the section, the portion of 
this outflow rate bound for destination r is given by 

vHt) 

Ti 

Note that the outflow rate given above depends only on 
the state of section i of the sub-network associated with 
destination r, yl(t). As a result, the dynamics for each 
destination are decoupled from all other destinations, 
which allows for a relatively straightforward extension 
of the solution of the single destination problem to the 
multiple destination problem. This is in contrast with 
the nonlinear outflow model, given in (2), which addi- 
tionally depends on the full state of section i , x,(t). 

The dynamics of the sub-network associated with desti- 
nation r can be written as 


Constraints will be given to ensure that conditions re- 
quired for Positivity , Conservation , and capacity con- 
straint satisfaction are satisfied at the end points of each 
interval of length At, i.e. at t ^ = kAt, k = 0, ..., K. We 
then give a method of generating continuous time rout- 
ing parameters which ensure that these conditions hold 
throughout the intervals. 


y r (t) = A r ( V r (t))y r (t) + B rT d flt), 

y r (o) = y r 0 , 

where 

A r {r] r (t)) = A 0 + ( 10 ) 

i&M r j&Ol Tl V ' 


Throughout the development of this solution method, we 
assume that the inflow rates d r s (t) are piecewise constant 
over the intervals of length At for all r £ [77] and s £ [5] . 
Inflow rates which conform to this assumption may be 
conservatively generated from more general inflow rates 
by taking d r s ( t ) as the maximum of the true inflow rate 
over the interval [tk,tk+ 1 ] where k = [^J. 


and A 0 




and 


B r {t) 


bsiit) bsnit) 


( 11 ) 


3. 1 Multiple Destination Linear Outflow Model 

3.1.1 Model 

Here, we extend the solution to Problem 1 presented 
in [1] to address the multiple destination problem with 
linear outflow rates. To clearly distinguish the linear and 
nonlinear models, which will be used in conjunction to 
prove results in subsequent sections, we will introduce 
new notation to describe the linear system. Sub-network 
state vectors y r , full network state vector y, and routing 


We denote the solution to (9) with initial state y r ( 0) = 
for all r £ [77] and inflow rates d r s ( t ) for all s £ [5] and 
r £ [77], under the particular choice of routing strategy 
V(t), as y I, (t). 

3.1.2 Basic Control Design Objectives 

Positivity of the newly constructed linear system (9) is 
defined in the same way as it is for the general nonlinear 
system (3) in Section 2.2. Note that linear system (9) 
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can be generated from the more general system (3) with 
rj ) replacing (3 and outflow rate 



which clearly satisfies the assumptions of Claim 1. Thus, 
we can conclude that Positivity and Conservation hold 
for system (9) if C(rj r (t)) hold for all r £ [R]. 

3.1.3 Control Design 

We now give sufficient conditions which can be used to 
generate a time- varying routing strategy rj(t) ensuring 
that the solution to system (9) remains below a contin- 
uous capacity bound c{t), that is y v (t) < c{t). We focus 
on each sub-network separately, giving sufficient condi- 
tions to generate routing strategy rj r {t) ensuring that 
the state of sub-network r remains below a continuous 
capacity bound c r {t), that is y r {t) < c r (t). We then de- 
fine the capacity bound for the full network as 

c 0) = J^c r (t)- 

r— 1 

These conditions are stated in the following proposition. 

Proposition 2 Given a fixed value ofr£ [i?], letrj r (t) 
be given such that rfj{t) > 0 for all i,j £ J\f r and t > 0 
and let t K > c r {t) be a differentiable vector-valued map 
such thatc r (t) > 0 for all t > 0. Then, ify r g < c r (0) and 

A r ( r/ r (t )) c r {t) + B rT d r (t) < c r (t) 

for all t > 0 the solution of system (9) satisfies y r (t ) < 
c r {t ) for all t> 0. 


since £i(t) = 0 and ( t ) > 0 implies that cf ( t ) — yf ( t ) > 
0. Thus, £(i) > 0 for all t > 0, and consequently y r (t) < 
c r (t) for all t > 0. □ 

Next, we derive constraints which can be imposed at the 
end points of each interval and a method of interpolation 
of rj(t) between these end points which ensure that the 
conditions needed to apply Proposition 2 hold through- 
out the interval. 

We make the choice to use a piecewise linear capacity 
bound, c r (t ), for each destination r £ [i?]. We parame- 
terize this function as 

r <=0, t=t 0 , 

c r (t) = < C r (t k -i) + Atm r (t%_i), t = tk>t 0 , (12) 

[ c r (t k ) + (t-t k )m r {t%), t^t k , 

where k = and m r (t) £ is constant over the 
intervals {t k ,t k+ 1 ). 

The choice of a piecewise linear capacity bound c r (f ) al- 
lows for a straightforward interpolation of r/(t) between 
the endpoints of time intervals, while preserving the 
properties needed to ensure that Positivity and Conser- 
vation hold and the performance objective is satisfied. 
The method of interpolation is given in the following 
theorem. 

Theorem 3 Given a fixed value ofr £ [i?], let the capac- 
ity bound vector c r (t) be given as in {12) for all t £ [0,if], 
and yg < Cq. For each k £ {0, . . . , K — 1}, if there exist 
rf (ffc) and rf (tfc+i) such that constraints C(r) r (t k )) and 
C{r) r {t k +i)) are satisfied and 


PROOF. Let £{t) = c r {t) — y r (t). We show that, for 
all i £ A f r 


at) = o 
m > o 


=► at) > o, 


which, since £(0) > 0, implies that £(t) > 0 for all t > 0. 
By definition of £ 


at) = m-m 
at) 


> - 


c r At) ® 


r-ieO 7 : 


'■? s =i 

s 


E i;.(‘)v-Ew) 

j.i&or s=i 


at) 

Ti 


E *?*(*) 


c r At)-y r At) 


i-ieo r . 


> 0 


A r (V r (tk)) c r (t k ) + B rT d r {t k ) < m r {t+) 
a r {tk + i)) c r (tk+ 1 ) + B rT d r (t k+ i) < m r {ti) 

for all i £ M r , then the parameters ri r (t) defined by 
Vij(t) = 

(! - ^r) + t -&ygt k+1 y i (tk+ 1) (14) 

(i - atk) + ^<(t fc+1 ) 

for all i,j £ M r and t £ [0, H\ where k = L^tJ; sat- 
isfy constraints C{rf {t)) thus ensuring that System (9) is 
Positive and Conservative for all t £ [0, H\. In addition, 
the resulting solution y r (t ) satisfies y r {t) < c r {t) for all 
t £ [0 ,H}. ' 

The proof of Theorem 3 is given in [1], The proof relies 
on the fact that rjf (t) is a convex combination of rfj{t k ) 
and rilj(tk+i) and y lj{t)c\{t) is a convex combination of 
r]ij{tk)ct{t k ) and ??t(tfc + i)ct(t fc+1 ), and thus constraints 
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C(r] r (t)) which are linear in 771 (f), and constraints (13) 
which are linear in are satisfied for all t £ 

[0 ,H\ whenever they are satisfied at the end points, t k , 
for ft = 0, . . . , K. 

Constraints (13) can be transformed into linear con- 
straints by introducing variable z\ 3 (t) to substitute for 
the nonlinear terms, for all i, j £ J\f r and all 

re [R], These constraints are given below, denoted by 
$ r : 


c r (t Q ) > y r 0 
for all ft e {0, . . . , K} 
c r {t k ) > 0, 

zl 3 (t k ) > 0, V*,j£ Af r , 
z r ij{tk) < cl(t k ), V ij e A A, 
zl 3 (t k ) = 0, V j e J\f r \Ol , V i e A/'' 
E”=i4» = <S(i*),VteAfVS£, 

" Ej:ieO 


Cj(tk) 


Z ^ + E S s =iKA(tt)<rnl(4), 


for all k e {0, . . . , K — 1} 

C r {t k + 1 ) = c r (t k ) + Atm r (t^), 

+ Ej:ieo? ^ ~ ~ + ZLi KA(t- k+ i) < ml(4) 


For each re [R], the values of r/ r at the end points of 
the intervals can be recovered as 


Vljitk) 


Zjj(tk ) 


(15) 


for all i, j e Af r , ft = {0, , K}. These values can then 
be interpolated between endpoints according to (14). 

A solution to Problem 1 can be generated by finding a 
feasible point of the linear constraints for all r £ [A] 
which also satisfies c (t) < c (t). If this problem is not fea- 
sible, we would like to give a recommendation on how to 
adjust capacity constraints c (t) in order to find a rout- 
ing solution satisfying these adjusted constraints which 
are, in some sense, close to the original desired con- 
straints. For this reason, we introduce piecewise constant 
adjusted constraint c(t). We then propose a solution to 
Problem 1 as finding a feasible point of 4) r for all r £ [A] 
which also satisfies c (t) < c(t) with the additional con- 
straint c(t) > c(t) while minimizing the integral of the 
difference between c(t) and c (t). When c(i) = c (t) the 
given capacity constraints are satisfied. 


problem: 


K - 1 N 

r mi. 11 - J2 ( 6 * (*fc ) - (4 )) At 

c r ,m r ,c,z z L ' 
k — 0 i= 1 

subject to $ r ',V r £ [i?] 

for all k £ {0, ... , K} (16) 

R 

E cr (^) - min {c(tfc), c(t+)}, 

r= 1 

c(t+) > c(t+). 

As defined earlier, the resulting capacity bound on the 
total amount of material in each section is 

c 0) = <?(*)■ 
r— 1 

Given that, under the routing strategy generated from 
a feasible solution of LP problem (16), the solution to 
system (9) satisfies y r (t ) < c r (t) for all t £ [0 ,H], it 
follows that y (t) < c(t) < c(t) for all t £ [0, H], 

Note that the objective function of LP problem (16) can 
be any convex function of the control design variables. 
For instance, a different positive weight could be used 
for each i to indicate that going over capacity in certain 
sections is worse than going over capacity in others. Al- 
ternatively, including the term 

X> r (^) 

r= 1 

in the objective function would provide an incentive for 
clearing material out of the network by the final time 
step. 

3.2 Multiple Destination Nonlinear Outflow Model 
3.2.1 Control Design 

Here, we leverage the control design method developed 
in Section 3.1.3 to solve Problem 1 with general nonlin- 
ear outflow rates. The decoupled dynamics of the linear 
system allowed us to treat each sub-network separately. 
However, in the nonlinear system, outflow rates of each 
section of a given sub-network depend not only on the 
associated state of the sub-network, but also on the full 
state of that section. Specifically, if we were to formu- 
late a proposition for the nonlinear system which paral- 
lels Proposition 2 for the linear system, the proof would 
require us to lower bound the term 


The problem of finding capacity constraints c and ca- 
pacity bound c can then be written as the following LP 


<$(*) 

Cj(t) 


A?( c j (0) 


Xj{t) 


Ml(xj(i))- 
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Unfortunately, we have no way to bound this term within 
the given framework. 

Rather than solving the nonlinear outflow rate problem 
directly, we take the solution method for the linear out- 
flow problem and add additional constraints on the rout- 
ing parameters which force recirculation within each sec- 
tion of the network. This recirculation effectively slows 
down the linear system so that the controlled outflow 
rate of each section of the linear system is at or below 
the corresponding uncontrolled nonlinear outflow rate. 
These conditions are formally stated in Theorem 4 be- 
low. 


Theorem 4 Let c r (t ) be defined as in (12) with m r (t £) 
for k = 0, . . . , K and c’ 0 > Xq given for all r £ [1?] . Let an 
outflow function y satisfying assumptions (1) be given. 

DefineTi = for alii £ [iV]. If,forallr £ [ii], 

there exist T) r (t k ) for k = 0, . . . , K , such that constraints 
C(r] r (t k )) are satisfied and the following constraints hold 
for k = 0, . . . K — 1 


A r (r, r (t k ))c r (t k ) + B rT d r (t+) < m r (t+), 
A r (r, r (t k+1 ))c r (t k+1 ) + B rT d r (t+) < m r (t+), 


(17) 


il r ii {t k )>ri i {k),\/i£U r 1 
rf ii {t k+ 1 )>r li (k),\/i£N r , ( j 

where A r {rf{t k )) is defined as in (10), B r is defined as in 
(11), and the fixed scalar value y^k) is chosen such that 


V i( fc ) > 


max 

tk |-i 


Tillj(Ci{t)) \ 

C i(t) ) 


for each i £ N and k £ [0, . . . , K— 1] , then the closed-loop 
system (3) under the decentralized time-varying state 
feedback control policy 


m) 

W) 




i(t) 




V i £ Af r 


ol 3 (t) X f Vz,je.AU\z^j 


(19) 

(20) 


withr) r (t) interpolated as in (14), has the following prop- 
erties 


(i) The solution (t) of system (3) is identical to the 
solution y* 7 (t) of system (9) with routing parameters 

v(t), 

(ii) Positivity holds for system (3), 

(Hi) Conservation holds for system (3), 

(iv) The solution x' 3 (t) of system (3), satisfies (t) < 
c(t) fort £ [0,77]. 


PROOF. We first note that the feedback control pol- 
icy defined by (19) and (20) was chosen so that the non- 


linear system (3) in closed-loop is identical to the linear 
system (9) with routing parameters rj(t). To see this, we 
substitute the expressions for (3(t) given in (19) and (20) 
into (3), which describes the dynamics of each section 
i £ J\f r of the sub-network associated with destination 

r e [R], 


W=-^ /*(*(*)) 


+ 


^(t)‘ 

1 - (! - ViM) 

E Um 


Xi(t) 


J 

s 


TiHfixfit)) 

Xj(t) 


*j(t) 


Xi(t) 

m 


7b(xi(t)) 


J Vj(Xj(t)) 


Hfixfit)) 


+ 


S— 1 

x r j(t) 

Xj(t) 


X-i [t ) Ti 

s 




S=1 


i-izo J 

. v- 


E itw 

j-ieo 


s 


which is the same as the description of the dynamics 
of each section i £ J\f r of sub-network r £ [Ii] given 
in (8). Thus, if y r (0) = 0^(0) = x r Q for all r £ [R] 
then x r (t) = y r (t) and (z) holds. Additionally, since 
inequalities (17) are the same as those in (13), and rj r (t) 
is interpolated between the endpoints t k of each interval 
by (14), for k = 0, . . . AT, Theorem 3 can be applied to 
conclude that Positivity holds for system (9), that is 
y r (t ) > 0, for all r £ [ii], and y(t) < c(t), for t £ [0, H\. 
Given (z), these properties also hold for system (3), thus 
(ii) and (iv) hold. 


Finally, we will show that Conservation holds by verify- 
ing that C((3 r (t)) hold for all r £ [ii]. From Theorem 3, 
we know that constraints C(rj r (t)) hold for all r £ [ii]. 
It is easily seen that C 2 (/3 r (t)) holds whenever C 2 (rj r (t)) 
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holds. To check C 3 (/3 r (£)) and Ci([3 r (t)) we compute 


E W)= E W)+Pii® 

jzov jeoi, 

= E + 1 


jeoi, 

3& 


- (i-^O)) 


Xi(t) 


= E »?«(*); 


Xi(t) 


jeOf 

Xj(t) 




+ 1 - 


iW 




E »?«(*) - 1 ) + L 

U'eor 


(21) 


Given that constraints C^{rf{t)) are satisfied, that is 

E ^(t) = i,vi£Ary>£, 

jeoj 

we can conclude from (21) that constraints Cs( y (3 r (t)) 
hold, that is 

E W = l,V*e-ATV£. 

ieor 

And given that constraints Ci(r) r (t )) are satisfied, that 
is 

E olj(t) < Vi £ 

jeo.r 

we can conclude from (21) that constraints C4(/3 r (£)) 
hold, that is 

E W) < ^ e «*£• 

ieor 


for all i £ A/’ 1 ’, r £ [i?], t £ [0, TV] implies that (22) holds. 
In order to show (23) we use the facts that fit is concave, 
/Lti(0) = 0, 0 < Xj(t) < Cj(t) which follows from (ii) and 

( iv ), and x 4 (f) = ^jc :*(£) + (l - x 0, therefore 

Mi(xi(t)) > ~^fj,i(ci(t)) +(l~ Mi(0) 

Mi(xi(t)) > ^|/i*(ci(t)) 

Xi(t) ~ a(t) 


From (14) and the fact that c[ (t) is piecewise linear, we 
have 


ViS)c r i(t) = ^1 - t -^^j Vii(tk)ci{t k ) 

+ ^^(< fc+ i K(t fc+ i) 

> 4 (%*(*) 

where the inequality follows from constraints (18). Thus, 
> v(k) where k = [^J. Given the definition of 
rj.(k), it follows that (24) and thus (22) hold and conse- 
quently that (4) holds for /3 V.(t) defined as in (19). 

Since constraints C(/3 r (t)) hold for all r £ [i?], we can 
conclude that Conservation holds for system (3) in 
closed-loop, thus (Hi) holds. □ 


The idea motivating the development of Theorem 4 is 
that the outflow rates of each section of linear system (9) 

with Ti = for all i £ [TV] can be restricted 

through constraints on recirculation so that the linear 
system has controlled outflow rate less than the corre- 
sponding uncontrolled nonlinear outflow rate, that is 


Clearly constraint Ci(/3 r (t)) holds when i ^ j for all 
r £ [R], In order to show that C\(j3 r (t)) holds when 
i = j, we must show that 


Vu(t) > 1 


TiHi(Xi(t)) 

x»(i) 


(22) 


for all i £ A f r , r £ [i?], t £ [0, H], First, we show that 


tii(xj(t)) Hi(cj(t)) 

Xi(t) ~ c i(t) 


(23) 


whenever Xj(t) < c i(t) and thus 


E (!-<(*))— <*&&)) (25) 

r= 1 

where the left hand side is the total controlled outflow 
rate of section i of the linear system. Constraints (18) 
ensure that (25) is satisfied. Given that r].(k) is chosen 
based on the upper bound c *(£) of the state y i(t) itself, 
these constraints are conservative. Recirculation beyond 
that required to satisfy constraints (18) may be used to 
satisfy the control design objective. Additional recircu- 
lation in the linear system beyond that required to sat- 
isfy (25) with equality is translated into recirculation in 
the controlled nonlinear system. 


V u(t) > 1 


TifJ,j(Cj(t)) 

c i(t) 


(24) 


We would like to write an LP problem, similar to LP 
problem (16), which can be used to generate parameters 
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which satisfy the assumptions of Theorem 4. Constraints 
(17) and C(r] r (tk)) for k = 0, . . . , K — 1 are ensured by 
constraints <I> 7 ’. The constraints 

4(M > V 

ZiSk+ 1) > r?.(fc)c-(t fc+ i) 


will ensure that constraints (18) hold. 

What remains is to find an appropriate choice of r/.(k) 
for all i £ [./V] and k = 0 , ,K. If a piecewise linear 
upper bound c(t) on c(f) were known, we could define the 
desired lower bound on recirculation parameter r/ 7 )- ( t ) as 


V S k ) 


max 

t k fe + 1 



Ti 


/bfe(s)) \ 
Cj(s) ) 


(26) 


where k = |_StJ- Given the fact that c i{t) is defined to 
be linear between t k and t k +i for k £ [0, . . . , if — 1], and 
thus 

c i(t) < max Ci(s), 

sG{tfc ,tk+ 1} 

we could then proceed as in (23) to conclude that the 
expression in the right hand side of (26) achieves its 
maximum value at an end point of the interval. Since c is 
a design variable, we do not have an a priori upper bound 
available. Instead, we explicitly introduce c, use it to 
define rj and add constraints to ensure that c(t) < c(t). 


The problem of finding capacity constraints c and ca- 
pacity bound c can then be written as the following LP 
problem: 


““ 51 51 (g(^) ~ )) At 

k — 0 

subject to 

for all k £ {0, . . . , K} 
c (t k ) < nhn{c(f"),c(t+)}, 
c (tt) > c(i+), 

R 

^2c r (t k ) = C (tk), 
r— 1 

C (t k ) < c {t k ), 

Ziiih) > %(k)c^(t k ), v i £ A r r y r£[R\. 


(27a) 


(27b) 

(27c) 

(27d) 

(27e) 

(27f) 


3.2.2 Iterative Solution Method 

Here, we present an iterative algorithm which can be 
used to solve Problem 1 by successively solving instances 
of LP problem (27). At the beginning of the algorithm, 
the upper bound c(t) on the capacity bound c(t) is set 
infinitely high, which forces full recirculation, that is 
Vi (tk) = 1 for all i £ [IV ] 7 r £ [R] and k £ [0, . . . , K — 1] . 
After LP problem (27) is solved, an appropriate finite 
value of c(i) can be chosen. At each subsequent itera- 
tion, c (t) is decreased and, correspondingly, required re- 
circulation is also decreased. 

The algorithm terminates when the maximum difference 
between c(t) and c (t) falls below some specified thresh- 
old. This termination condition ensures that the imposed 
recirculation constraints are reasonable with respect to 
the recirculation constraints required based on the re- 
sulting capacity bound c(t). 


Algorithm 1 Iterative refinement 


p — 0 

Ci(t fc ) = oo, for i £ [A], k = 0, . . . , K 
r/.(k) = 1, for i £ [N], k = 0, . . . , K 
Solve LP problem (27) 
n(p) {t k ) = c(tk), for k = 0, . . . , K 


9 

10 


6i P) (b=) = max (y max , c - p) (tfc)) , for i 6 [N], 
k = 0, . . . , K 

7: rtf\k) = max te{liA+i} (l - , 

for i = 1, . . . , IV, k = 0, . . . , K 
repeat 

p = p + 1 

c i P) (tfc) = c- p_1) (tfc) -7 (c! P-1) (ffc) -c^ _1) (tfc)) , 


11 : 

12 

13 

14 

15 

16 


for i = 1, . . . , N, k = 0, . . . , K 
rzP(fc)= max te{tfc tfc+i} ^ 
for i = 1, . . . , N, k = 0, . . . , K 
c (tk) = c (p) (t fc ), for k m 0, . . . , K 
r)(k) = ?/ p )(fc), for k = 0, . . . ,K 
Solve LP problem (27) 
c (p) (b=) = c (t k ), for k = 0, . . . , K 


1 ei(c«(t))\ 

1 Tl f n(t) ) 


until max 


i,k (c- p) (ifc) -cP(t fc )| 


< £ 


User defined values e > 0 and 0 < 7 < 1 affect the 
termination and speed of convergence of the algorithm, 
respectively. We define y max as the greatest amount of 
material that could occupy any section in the network by 


The specific choice of c (t) is important. If c (t) is chosen 
too high, then ryL(t) is forced to be high, resulting in a 
solution requiring more recirculation than necessary. If 
c (t) is too low, LP problem (27) may be infeasible. In 
Section 3.2.2 an algorithm is presented in which c (t) is 
adjusted iteratively allowing an arbitrarily tight upper 
bound on c(f) to be converged upon. 


y max = J2 l T y r 0 + 55 E E d r s (t+)At. 
r—1 r= 1 s=l k—0 


In the following proposition, we show that Algorithm 
1 is an effective method , that is, it terminates and is 
consistent. 
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Proposition 5 Given a network described as in ( 3 ) with 
section outflow rates p satisfying assumptions (1) and 
a piecewise constant desired capacity constraint c (t) for 
t G [0, H\, e > 0 and 0 < 7 < 1, the following properties 
hold 

(i) LP problem (27) is feasible when solved at Step 4 of 
Algorithm 1, 

(ii) LP problem (27) is feasible when solved at Step 14 
of Algorithm 1, 

(Hi) Algorithm 1 terminates. 


PROOF. First, we show that LP problem (27) is fea- 
sible in Step 4. At the first iteration 77 (&) = 1 for all 
i G [N] and k G [0, . . . , K — 1], This implies recircula- 
tion of all of the outflow of each section of the network. 
To construct a feasible solution, make the choice of: 


(27) at iteration p—1 and 7 < 1. Thus, any c,; feasible 
at iteration p—1 satisfies (27e) at iteration p. 

We now focus on constraint (27f) and parameter 77. 
Note that, for any feasible solution of LP problem (27), 
^{tk) < c- p-1 ^(ffc). Given 0 < 7 < 1, it follows from 
the definition of £( p_1 ) given in step 10 of Algorithm 
1 that c[ p \t k ) < c- p_1) (tfc). Also note that constraints 
<f> r ensure that c (t k ) = )Cf=i ° r (tk) > 0, and constraint 

(27e) ensures that £ (tk) > c (tk), thus c[ p \t k ) > 0. 
Given this, the fact that pi is concave, pi( 0) = 0 and 
that 

hM p) (t k )) > ^(c| p ~ 1} fa)) 
c i P) (t k ) c i P_1) (t fe ) 

whenever c[ p \t k ) < c- p_1 ^(ffc), as shown in the proof of 
Theorem 4, we have, 


Z ij 


(t k 


c r i(t k ),i = j, 
0, j 7^ i, 


m i(t£) = max \^b r si d r s (tl)^b r si d r s (t 


k+ 1) 


\s = 1 


«=1 


for all i,j G A f r , r G [i?], and k G [0, . . . , K — 1], Choose 
c r (t 0 ) > Xq for each r G [R] and feasible values for c r (tk) 
for all r G [R] and k G [0, . . . , K — 1] can be generated 
using (12) given m r (t k ). It is then straightforward to 
verify that constraints are satisfied for all r G [R]. 
Given that initially, £ (tk) = 00 for all fc G [0, . . . , K — 1] 
it is easily seen that constraints (27b) through (27f) also 
hold, and thus ( i ) holds. 


/ 1 AMi(c- P_1) fe)) 

- " #-%h) 

= p( p - 1) (fc). 

It follows that constraint (27f) can be satisfied at itera- 
tion p by z and c which are feasible at iteration p—1, 
recognizing that any feasible 2: and c are non-negative. 
Thus, we can conclude that any feasible solution of LP 
problem (27) at iteration p — 1, p > 0, is also a feasi- 
ble solution of the LP problem at iteration p and conse- 
quently (ii) holds. 


Given a feasible solution to LP problem (27) at iteration 
p — 1 for p > 0, we show that this solution is also feasible 
at iteration p. At each iteration of Algorithm 1, c and rj 
are the only constraint parameters of LP problem (27) 
which change. In order to show that LP problem (27) is 
feasible at each iteration, we show that the updated val- 
ues of c and 77 lead to a relaxation of constraints in the 
LP problem compared to the previous iteration. Param- 
eters c and 7j appear in only one constraint each, that 
is, constraint (27e) and constraint (27f), respectively. 

We begin by focusing on constraint (27e) and parameter 
c. For all i G [ N ] and p > 0, we have 

c [ P ~ 1] (tk) = c i P_1) (i fe ) - (cf _1) (t fe ) - <4 P_1) (f fe )) 

< c [ P ~ 1] (t k ) - 7 (c- P_1) (ifc) - c f -1) (t fc )) 

= c {p) (t k )- 

The inequality follows from the fact that cf~ l \tk) < 
ct 1] (t k ) given that c- p_1 ^ is feasible for LP problem 


We now show that Algorithm 1 is guaranteed to termi- 
nate. We focus on section i, time step tk . From a re- 
arrangement of the assignment of t r p \tk) in step 10, we 
have 


G (P 1] (tk) - c {p) (t k ) = 7 (c| p 1] (t k ) - c\ p 1} (t fe )) . 

(28) 

Constraints <f) r ensure that c- p ^(tk) < c- p ^(tk), thus 
the right hand side of the expression above is nonneg- 


ative and the sequence 



is non-increasing. 


This sequence is bounded below since constraints $ r en- 
sure that c (tk) = ° r (tk ) > 0, and constraint (27e) 
ensures that c(t k ) > c (tk)- Since sequence ^c j - P ' ) 


is monotone and bounded, it converges. Any convergent 
sequence is a Cauchy sequence and there thus exists an 
integer P ik such that for p > P ik 


|Si P 1] (tk) ~ Ci P) (t fe )| < 7 e - 


Substituting the expression for 


c - P 1) (tfc)-c- p) (tfe) given 
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iii (28) we have 

|Ci P_1) (ife)-c[ p_1) (tfc)| <e. 

It follows that for large enough values of p (i.e. p > 
max^fc Pik), the stopping criterion of step 16 will be sat- 
isfied and Algorithm 1 will terminate. □ 


Note that at each iteration of Algorithm 1, constraints 
are relaxed compared to the previous iteration. Thus, the 
cost of LP problem (27) decreases or remains constant 
from one iteration to the next. Physically, this means 
that the capacity constraint violation of the resulting 
closed-loop system decreases or remains constant at each 
iteration. 

Upon termination of Algorithm 1 (and, in fact, at any 
iteration of the algorithm), routing parameters r](tk) 
and 77 (f) can be recovered from (15) and (14), respec- 
tively. Under the decentralized time-varying state feed- 
back control policy defined by (19) and (20), system 
(3) satisfies the Positivity and Conservation constraints, 
and the solution x^(t) to system (3) satisfies x ,3 (t) < 
c(t) < c(t) for all f £ [0,17]. 

4 Application 

In order to illustrate the control design technique pre- 
sented in Section 3, this technique was applied to an air 
traffic control example problem. The network used in 
this problem is shown in Figure 1 presented in Section 
2.1. Network connectivity for the full network and sub- 
networks associated with destinations 1, 2 and 3 can be 
inferred from Figures 1(a), 1(b), 1(c), and 1(d), respec- 
tively. Note that, although not explicitly depicted, recir- 
culation is allowed in all sections of the network. 

Initial conditions and inflow rates are depicted graphi- 
cally in Figure 2. Inflow rates are constant for the dura- 
tion of the planning horizon of 3 hours. Flow rates and 
initial states are broken down by destination in Figures 
2(b), 2(c) and 2(d). That is, in Figure 2(b), the inflow 
rates specified are the inflow rates of flights with desti- 
nation 1, that is, &kdg(t) for each source, s = 1, 2, 3 and 
each section accepting inflow, i = 1,2, 3. The values in 
the boxes in figure Figure 2(b) indicate the initial num- 
ber of aircraft in each section with destination 1 , that is, 
mathematically, x)( 0 ) for all i £ A f 1 . 

The outflow rate, pi, as a function of the number of 
aircraft in the section is given in Figure 3. This outflow 
rate function is used for all sections of the network. For 
each section i £ [A] , 




(a) Full network. 



(b) Sub-network for destination 1. 



(c) Sub-network for destination 2. 



(d) Sub-network for destination 3. 


Fig. 2. Figures 2(a), 2(b), 2(c), and 2(d) indicate inflow rates 
and section initial conditions for all flights in the network 
and flights with destination 1, 2 and 3, respectively. The 
inflow at sources is indicated in the ovals on the left of each 
diagram. The value in each box represents the initial number 
of aircraft in that section. 



Fig. 3. Section outflow rate as a function of the total number 
of aircraft in the section. The same outflow rate function is 
used for each section of the network. 

Each section, except for section 14, has a constant ca- 
pacity constraint of 15 aircraft, that is, Cj(f) = 15 air- 
craft for all i £ [A]\{14} and all t > 0. Section 14 has 
the piecewise constant capacity constraint profile c.\^(t) 
pictured in Figure 4(a). 

Given these initial conditions and capacity constraints, 
Algorithm 1 was used to generate time-varying routing 
parameters to solve Problem 1. For this example, 7 = 0.7 
and e = 0.5. With these values, the stopping criterion 
was met after 16 iterations. At the final iteration, the 
cost associated with LP problem (27) was zero aircraft 
x hour, indicating that the given capacity constraint is 
not violated under the resulting routing solution. 

The given capacity constraint, and resulting capacity 
bounds and simulated state for section 14 are shown in 
Figure 4. The capacity constraint, Ci 4 (f), capacity bound 
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(a) Full state results for section 14. 



(b) Section 14 destination 1. 



(c) Section 14 destination 2. 



(d) Section 14 destination 3. 

Fig. 4. Figure 4(a) shows the capacity constraint 014 ( 1 ), ca- 
pacity bound C 14 (i) and state xi 4 (t) of section for 14 for the 
full network. Figures 4(b), 4(c), and 4(d) show the capacity 
bound cU(t) and state (f) of section for 14 for the sub- net- 
works associated with destinations r = 1, 2, 3, respectively. 



Fig. 5. Capacity constraint C 14 , adjusted capacity constraint 
C 14 , capacity bound C 14 and state X 14 of section 14. Note that 
Ci 4 is only plotted over intervals in which it differs from ci 4 . 

Ci 4 (i) = c i 4 (^) and full state, xi 4 (t), of section 14 

are shown in Figure 4(a). Capacity bounds c[ 4 (t) and 
state x\ 4 (f) are also given for r = 1, 2, 3 in Figures 4(b), 
4(c), and 4(d), respectively. 

In a second example problem, capacity constraints for 
section 14 were lowered in order to illustrate the capa- 
bility of adjusting the capacity constraints if the given 
capacity constraints cannot be satisfied. All other con- 
straints and parameters are identical to the above exam- 
ple. The given capacity constraint, and resulting capac- 
ity bounds and simulated state for section 14 are shown 
in Figure 5. 

Given these initial conditions and capacity constraints, 
Algorithm 1 was used to generate time-varying routing 
parameters to solve Problem 1. With these values, the 
stopping criterion was met after 15 iterations. At the fi- 
nal iteration, the integral of the difference between ad- 
justed constraint c (t) and given constraint c(f), calcu- 
lated as the cost of LP problem (27), is 1.28 aircraft x 
hour. Physically, this means that using this solution, the 
actual section count may exceed the constraint c(f) by 
no more than an average of 1.28 aircraft over a one hour 
time period. 


5 Conclusions and Future Work 

We addressed the problem of routing design for posi- 
tive compartmental systems with concave section out- 
flow rates to satisfy time- varying state constraints. The 
control design technique was presented as an algorithm 
which can be used to iteratively solve instances of an LP 
problem, successively tightening bounds on the closed- 
loop system state thus improving the quality of the so- 
lution. The resulting routing strategy is a time-varying 
state-feedback control strategy which guarantees that 
the state of the closed-loop system remains below the 
given capacity constraint, or an adjusted capacity con- 
straint if a solution to the given problem cannot be 
found. 
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Future work in this area involves the incorporation of 
uncertain capacity constraints. Such constraints arise in 
air traffic management problems in which section ca- 
pacity constraints may depend on weather conditions, 
which are not deterministic. A control technique which 
is robust to such uncertainty or reactive in real time to 
changing constraints would be advantageous for appli- 
cation in this field. 
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